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TUTORIAL 

Modeling and Simulation of Count Data 

EL Plani'2 

Count data, or number of events per time interval, are discrete data arising from repeated time to event observations. Their mean 
count, or piecewise constant event rate, can be evaluated by discrete probability distributions from the Poisson model family. 
Clinical trial data characterization often involves population count analysis. This tutorial presents the basics and diagnostics of 
count modeling and simulation in the context of pharmacometrics. Consideration is given to overdispersion, underdispersion, 
autocorrelation, and inhomogeneity. 
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OVERVIEW 

Pharnnaconnetrics,^ while having first been solely applied 
to continuous data due to its historical ties with pharnnaco- 
kinetics^ and its methodological complexities,^ now com- 
monly includes analysis of discrete type data.^ This tutorial, 
intended for pharmacometricians familiar with basic concepts 
in population modeling, simulation, and model-based drug 
development,^ aims at presenting the foundations of count 
data analyses and builds on a recent tutorial on time to event 
models.^ After defining count data and alternative analysis 
approaches, the main count models will be described with an 
emphasis on their assumptions, which will be completed by 
considerations in the context of drug development. 

COUNT DATA THEORY 
Count data definition 

Count data are generally defined as numbers of events per 
interval. For count data to be recorded, events need to be 
specified beforehand, and then recognized and counted dur- 
ing the experiment. The events should, in essence, be able to 
and, in practice, tend to occur several times in order for the 
gathered data to be qualified and identifiable as count. Math- 
ematically, counts are non-negative integers, since an event 
cannot happen an incomplete or a negative number of times. 
Naturally zero event/count is a possibility; certain entities 
may even display particularly many null observations, which 
phenomenon will be addressed later in this tutorial. Techni- 
cally, there is no upper boundary to a count, because there 
can theoretically be close to an infinite number of events 
taking place. 

Count data are not restricted to biomedicine; in fact, numer- 
ous statistical advances in the domain took place within other 
fields. Early applications ranged from the number of convic- 
tions pronounced by jury panels per year to the number of 
phone calls received at a call center per minute, through 
the number of cavalrymen killed by horse-kicks per year. In 
preclinical or clinical trials, counts are typically reported for 
each subject, e.g., number of epilepsy seizures per month, ^ 
number of incontinence episodes per week,^ or number of 
gastrointestinal symptoms per day.^ Such individual data are 



assumed in this tutorial, therefore finding their application 
within the context of population analysis.^° 

Underlying hazard process 

Since count data are intrinsically event frequency mea- 
sures, the inherent connection with repeated time to event 
is evident. Events are counted within time intervals essen- 
tially when their exact time of occurrence is ignored, usu- 
ally for convenience and practical reasons. For simplicity, 
both for data record and for data analysis, time intervals are 
generally — although not systematically — of fixed length, e.g., 
months, weeks, or days. The shorter the time intervals, the 
less summarized the occurrences are. As in survival analy- 
ses, events are presumed instantaneous and assumed to 
arise only one at a time, i.e., the probability that an individual 
experiences two events at the exact same instant is zero. 

In consonance with the analogous event data, count data 
are driven by an underlying hazard. Hazard is the instan- 
taneous rate of the events.^ Determining whether and how 
this hazard varies with covariates, including treatment effect, 
is typically the aim of the data analysis. When time affects 
the hazard, i.e., the hazard is not constant, summarizing the 
information in counts per interval may dilute the effect. The 
extent of this dilution depends on the pattern of the time effect 
and its coincidence with the time interval limits. Nevertheless, 
in studies where the goal is to most accurately character- 
ize a time-varying hazard, a repeated time to event analysis, 
although slower, will be more appropriate than a count analy- 
sis. The exceptions to this principle are if the events are rare 
enough and if the time intervals are short enough for only one 
or zero count to be recorded. 

Noncount data 

Some types of data are sometimes unsuitably identified as 
event count data, such as grouped binary data or ordinal 
data. Grouped binary data arise from counting the number 
of successes or failures resulting from a certain number of 
binary outcome processes, i.e., Bernoulli trials. A binomial 
model, handling the number of tests performed and there- 
fore the implicit maximum, appears more adapted to analyze 
grouped binary data. A Poisson model can nevertheless 
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be appropriate to describe binomial data if the number of 
experiments is large, since the two models are mathemati- 
cally equivalent for an infinite number of repetitions. Regard- 
ing ordinal data, ordered categorical models are the suitable 
models for the evaluation. 

On the other hand, true event count data are often han- 
dled with models for a different data type. For instance, they 
may be categorized, i.e., dichotomized (event, no event) or 
made ordinal (low, medium, and high frequency of events). 
This practice, intended to simplify the process, risks informa- 
tion and power loss. Another alternative approach applied to 
count data is to treat them as continuous, which does not 
involve data transformation. A continuous model, despite not 
addressing the integer nature of the data, may be successfully 
applied to count data, predominantly when the distribution is 
wide with high counts. These characteristics commonly apply 
to events in a volume or an area, e.g., number of neutrophils 
per microliter of blood or number of bacteria per square centi- 
meter of a Petri dish. It should be noted that the classic count 
methods can also be used for these types of observations 
nonetheless. 

Count data visualization 

There are several — standard or not — ways to visualize count 
data, and a representative sample will be given in this tuto- 
rial. The familiar representation of the dependent variable — 
here the counts — vs. time (Figure 1 , left panel) is inevitable, 
although it does not carry the notion of the integer nature of 
the observations. This type of graph will, as for other data, 
inform about the variability between subjects as well as a 
potential time course. 

A histogram of the frequency of the counts (Figure 1 , right 
panel) is highly recommended, because of the indication it 
gives about the shape of the distribution. This is referred to as 
probability mass function (pmf) owing to the discrete nature 
of the probability distribution of counts, as opposed to the 
probability density function associated with continuous data. 



In Figure 1 , the histogram was intentionally rotated to exhibit 
how the two y-axes coincide, which is, however, not common 
practice. It can be highlighted that most variables included in 
the count model structure in fine govern the pmf. 

COUNT IVIODEL STRUCTURE 
The Poisson model 

The most important count model is the Poisson model. It is 
also the model to which all other count models presented 
in this tutorial converge, when their extra parameter goes to 
zero or infinity. It is named after the French mathematician 
who reintroduced the concept^^ a century after de Moivre 
published it.^^ 

The model describes the observations using the Poisson 
distribution. Accordingly, the random variable in the model is 
Poisson distributed. The term Poisson model used thereafter 
corresponds to this definition. 

An observation /during an interval yfor an individual / is 
said to follow a Poisson distribution when the probability that 
it takes values n = 0, 1 ,2, . . . can be expressed as in Eq. 1 . 

P{Yij = n) = ^-e-'' (1) 
^ ^ ^ n\ 

where X. is a parameter and ! the factorial function. 
Mean count parameter 

The parameter A (lambda) is the single parameter of the 
Poisson model in its simplest form. The expectation for the 
Poisson distribution E (V) is equal to A, and the parameter 
will therefore also correspond to the arithmetic mean of the 
counts occurring during a certain time. A is by definition a 
positive real number, i.e., it is not necessarily an integer. 

In pharmacometrics, parameters are often expressed as 
a mixed effect — a combination of a fixed (0) and a random 
(r|) effect. A common parameterization preventing negative 
values is given in Eq. 2. The introduction of interindividual 




Figure 1 Representation of count data. The left panel displays the time course of the counts in a 25-individual population, whereas the right 
panel reveals the probability mass function of the counts in the same population. 
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Figure 2 Properties of the parameter A. The plots represent how varying A values between individuals or groups of individuals result in shifted 
probability mass functions (left panel, with three distributions) and shifted mean counts (right panel, stratified on treatment). Varying A within 
individuals can also be revealed in the right panel with fluctuating means overtime, here represented with their 95% prediction intervals. 



variability (II V) — the randonn effect — results in A to be inter- 
pretable as the nnean nunnber of events per tinne interval for 
individual /. 



(2) 



As opposed to, for exannple, a nornnal distribution, the vari- 
ance of the Poisson distribution is not governed by an addi- 
tional paranneter, but by A. as well. As a result, an individual's 
A. value defines both the nnean and the variance of the dis- 
crete probability distribution for the observations of this spe- 
cific individual, as in Eq. 3. 

A,=E(\^) = Var(\</) (3) 

This feature of the Poisson distribution is illustrated in 
Figure 2 (left panel) with three pnnfs corresponding to three A 
values. It can be observed that as A increases: (i) the nnean of 
the distribution shifts and the probability of zeros decreases, 
(ii) the function is less skewed and approxinnates a nornnal 
distribution, and (iii) the variance increases and higher counts 
are encountered. 

Poisson model properties 

The equality between nnean and variance of the counts 
experienced by an individual constitutes an innportant prop- 
erty of the nnodel through its distribution. Manned equidis- 
persion, it should be investigated in the data exploration 
and disclosed in the nnodel report. Though restrictive on the 
statistical level, this trait naturally applies to nnany in vivo 
phenonnena. 

Another characteristic intrinsic to the Poisson nnodel is 
independence between events or nnennoryless process. 
Connnnon for discrete probability functions, it nneans that 
the occurrence of an event does not alter the probability of 



another. While adnnissible in idealistic cases such as rolling a 
die and counting the nunnber of occurrences a 6 is obtained, 
this assunnption nnay generally not be fulfilled when dealing 
with bionnedical synnptonns. 

The third Poisson property, presented here and discussed 
elsewhere in the tutorial, is the Constance of the event rate 
within tinne intervals, also known as honnogeneous process. 
On account of an evident lack of infornnation at this level, this 
aspect is inherent to the data collection. While echoing the 
consequences of low-resolution tinne intervals, the feature 
can be perfectly acceptable in nnany cases. 

Exponential model family 

Considering A as a piecewise constant rate further enables to 
establish the link between the Poisson nnodel and the hazard 
nnodel. 

The sinnilarity between the two exponential nnodels can be 
depicted, for instance, by connputing the probability of zero 
event to happen within an interval according to the Poisson 
nnodel (Eq. 4) and the sanne probability given by the tinne-to- 
event approach when the hazard is known to be constant 
(Eq. 5). 



P{Yij = 0)- 



n! 



-A/ 



0! 



-A/ 



■- e 



-A/ 



P(Yij = 0) = S(tj) = e 



) du 



(4) 



(5) 



The analytic solution of the integral in the case of a constant 
hazard nnakes nnore apparent that A, the variable describing 
the nunnber of events within an interval in Eq. 4, corresponds 
to the product in Eq. 5 of the rate of the events and the length 
of the interval. 
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The proof of the mathematical equivalence between a 
piecewise proportional hazards model and a Poisson regres- 
sion model was produced, along with demonstration of identi- 
cal estimates, independently by two teams in 1980J^'^^ 

The advantages of acknowledging the direct relationship 
between A and a constant hazard reside in the meaning of 
the methods, the understanding of the results, and the inter- 
pretation of the conclusions. Stating the existence of a clear 
correspondence (Eqs. 4 and 5) between the two models 
also suggests shared robustness and accuracy in terms of 
inferences. 

Mean count function 

If A, in the standard Poisson model does not vary within time 
intervals, it can nevertheless vary between them. Since A 
represents the mean counts, it can be susceptible to time 
changes and is the parameter on which a potential drug effect 
will primarily be explored. These elements place A at the 
heart of disease progression and pharmacodynamics char- 
acterizations, which include in silico description of the time 
course (natural history or placebo effect) of the disease status 
assessment and the action of the drug (if achieving a measur- 
able response). The disease status assessment is in this case 
the counted events chosen as end point during, e.g., a clinical 
trial. The assessment is expected to undergo a response qual- 
ified either as a desired effect or as an adverse event, both of 
which can take the form of a stimulation or an inhibition. 

After establishment of the count structural model, model 
building can be carried out on A the same way it is for a con- 
tinuous model, in terms of a variability model and a covari- 
ate model. In pharmacometric models, A is typically a function 
that can be affected by time (/(f)), drug effect (g(PK.(t^), and 
covariate effects (a^X^,a^X^,...a^^). Eq. 2 may therefore be 
extended to, e.g., Eq. 6. 

^i = ee^^O-f(tj))0-g{PK^(tj))- 

The variable PK.(t^ stands for a pertinent pharmacokinetic 
model-derived exposure metric^^ driving the response. It is 
apparent that the exposure variable will account for the entire 
time interval and, therefore, is likely to be a summary metric 
(e.g., steady-state concentration or average area-under-the- 
curve) as opposed to momentary concentration. 

Note that setting the functions f and g (in Eq. 6) to asymp- 
totically approach zero allows the mean count to tend to zero 
while keeping it from becoming negative. Prevalent expo- 
nential decay and sigmoid E^^^ functions, respectively, can 
be successfully used in this context. To the same end, oth- 
ers may prefer a parameterization on the natural logarithmic 
scale such as, e.g., Eq. 7. The fixed effect of the intercept 0 
in this case will not be directly comparable with the one from 
the previous equation. 

\og(Xi) = 0 + iii + f(tj) + g(PKi(tj)) 

+ a^X^ + a2X2 + ... + aKXK ^ ^ 

Figure 2 (right panel) illustrates such a case with the mean 
response in a hypothetical population displaying a decline 
over the month, which is more pronounced in one of the 
two treatment groups. The ordinate axis clearly depicts A as 



no longer only a parameter but a multidimensional function 
varying with time and exposure and furthermore, between 
individuals. 

COUNT IVIODEL COIVIPONENTS FOR DATA VIOLATING 
POISSON ASSUMPTIONS 
Nonequidispersed events 

Equidispersion can be verified numerically by computing and 
comparing the mean and the variance of the observations of 
each individual represented in the raw data. A statistical test 
consists in calculating the distributed ratio between the 
variance and the mean, i.e., the index of dispersion, ^^'^^ and 
testing equality to one. However, in a context of time-varying 
mean count, only limited conclusions can be drawn from this 
practice. Therefore, such calculations should only be inter- 
preted along with other diagnostics and must not preclude 
further investigations. 

Graphically, when individual means and variances are 
plotted against each other, a linear relationship suggests that 
the condition is respected. Violations of the equidispersion 
assumption can be of two kinds: overdispersion or underdis- 
persion. Overdispersion corresponds to cases of a variance 
greater than the mean for most of the population, and under- 
dispersion to the opposite, a variance lower than the mean. In 
the mean-variance plot (Figure 3, left panel), these types of 
profiles are respectively situated above and below the iden- 
tity line. Again, this may be impacted by a nonconstant A, and 
thus, a stratification on time (e.g., by week in this daily record 
example) should be considered. 

Another diagnostic designed to indicate whether the 
observed distribution approximately compares with a Poisson 
distribution is the Poissonness plot.^^'^^ It is based on a quan- 
tity called count metameter, calculated from the frequency of 
each encountered count (Q, the factorial of the count values 
(n!), and the sum of all frequencies (A/). According to Eq. 8, if 
the data follow a Poisson distribution, the metric is linear with 
regard to the possible counts n, with an intercept -A and a 
slope log (A). 

fn = ^- e-'-' ■N^\og(f,) = n- log (A, ) - A, 

^ (8) 
- log (n !) + log (A/) ^ log [^^] = + ^ • log (A, ) 

This diagnostic (Figure 3, right panel), similar in concept to 
the Q-Q plot, however, seems to be more adapted to detect 
overdispersed distributions — more points not falling on a 
straight line — than underdispersed ones, the double Poisson 
in this example. The diagnostic power of the available plots is 
limited due to the stochastic nature of the count distribution 
as well as the potentially complex A profile. Therefore, one 
should not completely rely upon them and investigate alter- 
native models described in this section. 

Overdispersion or underdispersion circumstances may not 
always be explainable. The first case, denoting that the within- 
individual variability is especially large compared with the 
mean, i.e., that the same patient can experience few events 
during most time intervals and many at some other times in 
the same study, can be encountered in types of diseases pre- 
senting sudden outbreaks of symptoms. The second case, 
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Figure 3 Dispersion plots. The left panel is a mean-variance plot and the right panel is a Poissonness plot. Each dot represents metrics 
of all counts for one individual on the left, and a metric for each possible count of one individual on the right. Each line consists in a linear 
regression of the dots belonging to the same model on the left and to the same individual on the right, and the additional dashed line on 
the left panel displays the identity. The dots on the right panel also have their 95% confidence interval represented and are colored black 
if the linear regression did not include this interval. The plots were produced from simulations with different models: Poisson, negative 
binomial (NB), generalized Poisson with a positive dispersion parameter (GP (5 >0)), and double Poisson (DP). 



referring to notably low within-individual variability vs. nnean, 
i.e., honnogeneous records, nnaybe displayed by rather con- 
stant experiences and robust synnptonns, like in diseases 
qualified as slowly progressive. The cause of the dispersion 
phenonnenon nnay be identifiable, e.g., epilepsy counts nnay 
substantially differ between nnorning and evening and the 
binnodal distribution nnay be nnisinterpreted as overdispersed, 
or recurring nnigraines nnay be treated with rescue nnedication 
leading the count distribution to appear underdispersed. But 
in nnany cases, driving forces for the noise in the process are 
nnissing or unknown, and resulting overdispersion or under- 
dispersion seen to occur are left to be characterized in the 
absence of explanatory factors. 

Analyzing data not confornning to the equidispersion 
assunnption related to the Poisson distribution nnay require 
the use of candidate nnodels different fronn the Poisson 
nnodel. Several nnodels have to this end been proposed and a 
few are presented here. 

Negative binomial model. The negative binonnial (NB) 
nnodel,22'22 also sonnetinnes referred to as the inverse bino- 
nnial nnodel, is the nnain alternative to the Poisson nnodel 
in case of overdispersion (Figure 4, right panel). It uses a 
supplennentary paranneter, O (onnicron), responsible for an 
increased variability connpared with the nnean counts. This 
overdispersion paranneter belongs to the range O e (0,oo), 
but values below 1 are generally sufficient. The discrete prob- 
ability distribution (Eq. 9) still has a nnean of E(K) = A, while 
the variance is now Var(K) = A.(1+ O.-X). 



Although the distribution is not defined for C> = 0, it tends 
toward a Poisson distribution when approaching it. This 
nneans that by allowing the paranneter to vary between indi- 
viduals — by including a randonn effect — both NB-like and 
Poisson-like profiles can be captured. The NB nnodel was 
found adequate when applied to epilepsy data in several 
analyses^'^'^^ and to infectious diseases.^^ 

Generalized Poisson model. The generalized Poisson 
nnodeP'^^ is another nnodel used to account for deviation 
fronn equidispersion in data to be analyzed. The additional 
paranneter, 5 (delta), governs the dispersion of the pnnf. The 
distribution can result in an overdispersed pnnf if d is posi- 
tive (d > 0) and in an underdispersed one if 5 is negative 
(d < 0) (Figure 4). The linnits within which d exists are d e, 
nnax (-1,-A/ / A),1), where A is the largest positive integer 
greater than or equal to 4 that satisfies the inequality X+Ad> 
0. The paranneterization in Eq. 10 pernnits the nnean count to 
be equal to E(K) = A. and the variance to Var(K) = 

Clearly, the generalized Poisson nnodel reduces to the Poisson 
nnodel and therefore adjusts for equidispersion when 5=0. 
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Figure 4 Nonequidispersed distributions. The same distribution density issued from a Poisson model (A = 3) is plotted on both panels, together 
with two models for underdispersion (narrower distributions) on the left panel and two models displaying overdispersion (wider distributions) on 
the right panel. The explored models are double Poisson (DP) (v = 3), generalized Poisson (GP) (5= -0.75 and 0.50), and negative binomial 
(NB) (0 = 0.75). 



The flexibility of the nnodel to acconnnnodate the twin proper- 
ties of over- and underdispersion nnakes it an attractive candi- 
date. Nonetheless, the dependence between the lower bound 
of S and A, two paranneters to estinnate, introduces an incon- 
venience. A step-wise procedure can be used to circunnvent 
the difficulty by starting with 5=0 and then with the obtained 
fixed effect estinnate A setting nnax (-1, -A / 4) < (5 < 1. 

Double Poisson model. The double Poisson nnodeP^'^° 
is generally pronnoted in the context of underdispersion 
(Figure 4, left panel), since in contrast to the previous nnod- 
els, the higher the extra paranneter, the nnore underdispersed 
the distribution. Both fornns of dispersion can be nnanaged 
though, through the paranneter v (upsilon). The real nunnber v 
yields underdispersion when above 1 (v> ^) and overdisper- 
sion when below 1 (v<^). The nnodel gains in flexibility when 
IIV is associated with v. The discrete probability distribution 
generated by Eq. 1 1 has a nnean of approxinnately E(V^) = A, 
and a variance of approxinnately Var(y^) A, / v'. 



P(Yij=n) 



^-vi 



n\ 



1 2Vj/ij y VjXj 



(11) 



As its nanne suggests, the nnodel is based on the double 
exponential fannily. This extension of the Poisson nnodel is 
obtained as an exponential connbination of the Poisson den- 
sity with nnean A and of the Poisson with nnean n. In the case 
of D = 1 , the double Poisson nnodel collapses to the Pois- 
son nnodel. Technically, the double Poisson function only 
approxinnates a pnnf. The last elennent of the equation is a 



nornnalizing constant specified to nnake the sunn of the prob- 
abilities closer to unity. 

Other nonequidispersion models. Several other nannes, nota- 
tions, or paranneterizations exist for the nnodels described 
above that do not correspond to different distributions. The 
nnentioned nnodels were used to sinnulate density profiles to 
connpare underdispersed (Figure 4, left panel) and overdis- 
persed (Figure 4, right panel) data to a unique Poisson — 
equidispersed — set of data. While the nneans of all pnnf were 
chosen to be the sanne, the dispersion paranneters were set to 
illustrate the feature, which could be enhanced or dinninished. 

Other nnodels were proposed to handle overdispersion, 
such as nnixture nnodels. This alternative is based on the sus- 
picion of two (or nnore) underlying processes leading to an 
event. A nnixture nnodel can be developed with two (or nnore) 
X values, where the observed sunn of events has a larger vari- 
ability than either of the two distributions. 

Another frequent case of overdispersion is when it is due 
to an excess of a certain value. Zero counts are connnnonly 
observed in higher frequency in the data connpared with 
the annount predicted by the nnodel. A possible situation is 
a threshold needed to be crossed for positive counts to be 
observed, another being a poor connpliance when adverse 
events are counted. The Hurdle nnodel, built as a nnixture of 
two different densities, nnay rennedy to excess of zeros. The 
zero-inflated Poisson nnodel is a straight-forward substitute, 
where the two nnixture connponents, (i) zero and (ii) a Pois- 
son pnnf, are weighted by a factor corresponding de facto to 
the probability of not observing any event. The zero-inflated 
Poisson nnodel can be connbined with distribution functions 
other than the Poisson nnodel, e.g., a zero-inflated NB nnodel 
was explored in the analysis of antibody titers^^ to account for 
patients not developing innnnunity. 
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Figure 5 Serial correlated counts. A population displaying correlation between consecutive observations (Markov Poisson) is represented 
in comparison with a simulation from a Poisson model (assuming independence between counts). The left panel investigates the correlation 
between neighboring counts with the dots representing the counts and the solid lines the linear regressions (dash line is identity line). In 
the right panel four individuals are exhibited, with slowly varying profiles in the Markov population and highly fluctuating curves in the vanilla 
Poisson simulations. 



More nonnnixture nnodels presenting dual properties of 
over- and underdispersion are: the connpound generalized 
Poisson distribution, 22 the Katz System yielding the general- 
ized event count model, ^^'^^ or the Conway-Maxwell Poisson 
model. ^^'^^ 

A further reason for observing underdispersion is data 
truncation. Right truncation is observed when events above 
a fixed value {c) are not counted, whereas left truncation 
occurs when the censoring affects counts below a set value 
(c). The most common case of the latter is zero truncation. 
Models to be used in case of data truncation can be any 
of the above, yet their distribution must be truncated while 
retaining probabilities adding up to 1. Density truncation can 
be applied though division by the cumulative density function 
(cdf) within the range of possible counts, i.e., sum of proba- 
bilities Xn=o " ^) for right truncation or 1 - XI=c " 
for left truncation. This technique was also used 'to model 
score data.2^ 

Nonindependence between events 

Assuming that events happen in an independent manner 
may not be appropriate to describe disease symptoms. 
In epilepsy, the notion that seizures beget seizures is 
debated, but implies that seizures themselves are epilep- 
togenic. Besides events directly triggering more events, 
events clustered in high-intensity periods or low-intensity 
periods can denote an underlying — indirect — dependence. 
Migraines are believed to follow such an unknown process 
leading low-activity epochs to intercalate with peaks. Mictu- 
rition contains a psychology factor leading events to be arti- 
ficially more regularly spaced compared with if truly driven 
by a random process. 



Dependence between events can be diagnosed by scrutiniz- 
ing the data. A simple way to visualize dependence is to plot 
consecutive observations against each other (Figure 5, left 
panel). A linear regression line can indicate in which fashion 
their dependence occurs. Individual profiles (Figure 5, right 
panel) may also support independence invalidation, since 
within-individual dependence in biomedical observations is 
likely to be expressed as trains of similar values. Correlations 
may be involved if the displayed time courses are smoother 
than when generated under an independence assumption. 

In statistics, this property is named serial correlation and 
can be addressed in different ways. In pharmacometrics, 
approaches implemented to reflect dependence between 
observations include the first-order autoregressive (AR(1)) 
model, 22 correlating residual errors, and stochastic differen- 
tial equations,29 decomposing residual errors into system and 
measurement noise. However, both approaches have been 
described so far only for continuous processes. 

For discrete models, published methods for serial correla- 
tion mainly relate to Markov models. Observations in Markov 
models are assumed dependent on previous states. The previ- 
ous state in question is often the exact preceding one, referred 
to as first-order Markov model, but can instead or in addition 
implicate the observation before the last (i.e., second-order 
Markov), etc. Because the chain of states is discrete, i.e., time 
intervals are defined as a set of measurement occasions, the 
full terminology is discrete time Markov chain. 

Note that the preceding observation is known and can 
therefore directly enter the model. This is a concept differ- 
ent from that of hidden Markov models, where the state 
influencing the current observation is unknown. In a hidden 
Markov model, the latent system transitions between hidden 
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States, with the state where it is positioned at a certain time 
altering the probabilities for the observations at that time. 
This methodology was implemented for mixed-effects count 
models/^ 

Markov models are integrated as components to probabil- 
ity models through functions of a prior observation affecting 
the probability of a current observation (Eq. 12). In count 
models, most implementations are specified according to 
models considered by Zeger and Qaqish.'^^ Markov compo- 
nents contribute to modify the pmf mainly via addition to A, 
or through inflation-deflation probabilities (similar to the zero- 
inflated Poisson model). 

P(Yij = n\Yij_,) = f(n;XrXM) (12) 

Markov element predictors can be of any type: binary, cat- 
egorical, count, or continuous, since they can be the actual 
previous observation or a derived value. In pharmacodynamic 
count models, the two most encountered parameterizations 
are a continuous function of the previous count and a mixture 
function of the previous status no count/counts. Examples of 
the first case are a proportional E^^^ function of the preced- 
ing epilepsy seizure count^^ or a linear function of the past 
number of multiple sclerosis contrast enhancing lesions"^^ 
(or a multiplicative power function of hypothetical events in 
Figure 5). Still applied to epilepsy, results were presented 
where the last observation information was dichotomized 
depending on whether the seizure was present or absent.^^'^^ 

Nonconstant event rate within time intervals 

As expressed before, the concept of hazard^^ is key in this 
context, since this rate determines when the events, ultimately 
counted, occur. Applying a model from the Poisson family 
assumes that this rate is constant within the intervals where 
the counting is made. These stationary increments may not 
be suitable for all problems at hand. Time variations of event 
occurrence are likely to take place within all time intervals, 
but their capture may be more possible and crucial in some 
circumstances than others. 

In the case of safety studies, when the count variable is 
a rare adverse event, this time variation is seldom of impor- 
tance and often not feasible to capture, given the sparse- 
ness of the available data. Similarly, the description of cyclic 
fluctuations, like a circadian rhythm, where the length of the 
period coincides with the length of the interval, is likely not 
possible. However, prior or simultaneous characterization of 
a time-varying biomarker or a pharmacokinetic profile, sus- 
pected to influence the hazard, may be the ground for requir- 
ing nonconstant event rate within intervals. 

Relaxing the stationarity assumption corresponds to con- 
verting the homogeneous Poisson process into an inhomo- 
geneous (or nonhomogeneous) Poisson process. Ample 
research was performed in the branch of recurring events rate 
analysis. Inhomogeneous Poisson models couple (i) models 
with time-varying count intensity^^ to (ii) Poisson (or other) 
probability distribution models. In this case, the expectation 
of the sum of counts during a finite interval (Eq. 13) can be 
obtained from the sum of the expectations of each event. 

^Y^U-xj)) = \l^^ii^)du (13) 



Where X.(u) can be expressed as in Eq. 6 with the P/C(f) being 
for example the drug concentration (not averaged over the 
time interval any more) or a covariate X being time-varying. 

Ordinary differential equations are, therefore, required 
most of the time in this family of models. Nevertheless, ana- 
lytical solutions may be employed in some cases (e.g., the 
Weibull distribution). The parameter/function X acts as haz- 
ard which accumulates over time. The likelihood of positive 
counts is set to the chosen pmf, while the likelihood of zero 
count is equal to the survivor function. The hazard, cumula- 
tive hazard, and survival were represented in Figure 6. 

The literature reporting pharmacometric investigations 
using inhomogeneous Poisson models is relatively sparse. 
One of the studies is a sensitivity analysis made with regard 
to the data simplification to compare homogeneous models 
with an inhomogeneous model. In other respects, a novel 
approach was proposed connecting a third type of model: an 
ordered categorical model, to estimate the maximal severity 
achieved from nonconstant repeated categorical events sum- 
marized per time intervals."^^ 

One may note that the inhomogeneous approach relieves 
the considerations with regard to the regularity of the interval 
length, and an integration is done over elapsed time after 
the beginning of the study. Intervals of different lengths can 
nonetheless be efficiently handled in a homogeneous model. 
As indicated by Eq. 5, the interval duration can be taken into 
account such as Eq. 14. Hence, if symptoms are, for exam- 
ple, counted every day during stays at the clinic and every 
week when the patient is at home, the data can be analyzed 
without having to add daily counts into weekly ones. 

^ij=^i-{tj-tj_,) (14) 



COUNT IVIODELING PRACTICES 
Model fit software 

Model fit of count data can be performed in all major nonlin- 
ear mixed-effects modeling software, i.e., the pharmacomet- 
ric software NONMEM, Monolix, Phoenix, and S-ADAPT, and 
the statistical software WinBUGS, R, SAS, and Matlab. 

The language (e.g., NM-TRAN) used to code the probabil- 
ity functions naturally depends on the software of choice, with 
the Poisson (and oftentimes some other) distribution(s) being 
predefined in a few of the programs (e.g., WinBUGS, where 
one can write "dpois"). 

Furthermore, the factorial function (•!), which is part of 
all count models, is not systematically embedded in all 
programs. Approximations such as a refinement of the 
Stirling's formula or the one derived by Ramanujan are a 
solution to this. The gamma function (T-) can be handled 
the same way since it is directly related to the factorial func- 
tion. For most models though, the approximation of the fac- 
torial does not affect parameter estimation as it is constant 
in the likelihood. 

Another potential difficulty taking place during count mod- 
eling is numerical issues while dealing with high counts. This 
risk can be prevented and model stability increased by log- 
transforming the model or its components, e.g., the gamma 
function of the NB model. 
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Figure 6 Inhomogeneous model variables. The plot displays variables for two individuals. The "count" variable corresponds to observed data 
and the "hazard" to the function A. The "cumulative hazard" and the "survival" function are quantities involved in the estimation, and scaled 
here to share the same y-axis. 



Estimation method choice 

Maxinnunn likelihood estinnation is the prevalent paranneter 
estinnation nnethod in pharnnaconnetrics. Since count nnodels 
consist in probabilities, they are easily interpretable as sets of 
likelihoods, in the sanne way the other types of discrete nnod- 
els are."^^ Estinnation in population count nnodels is, therefore, 
connpleted through likelihood-based approaches, requiring 
approxinnation of the integral over the randonn effects. As in 
other cases, estinnating paranneters through nnininnization of 
-2 tinnes the log-likelihood presents advantages over nnaxi- 
nnization of the likelihood. 

Maxinnunn likelihood estinnation algorithnns nnay not perfornn 
equally well due to the challenging nonlinearity of the nnod- 
els. Accordingly, the use of the prinnary NONMEM estinna- 
tion algorithnn FOCE is prohibited and the alternative Laplace 
and sannpling-based nnethods reconnnnended. The evalu- 
Qtjon5o,5i Qf ^i^g perfornnance — accuracy and precision — of 
the algorithnns Laplace in NONMEM, stochastic approxinna- 
tion of expectation nnaxinnization in Monolix, and Gaussian 
quadrature in SAS indicated that these estinnation nnethods 
are appropriate for count analyses. Bayesian nnethods are 
also successfully used to estinnate count nnodel paranneters 
as seen on available training nnaterial.^^ 

IVIodel development workflow 

After appropriate graphical exploration of the data, knowl- 
edge about the data structure, profile, and trends is undeni- 
ably gained. This infornnation can highlight features to reflect 
in the nnodel structure as well as to use in guiding nnodel 
building. 

The nnodel developnnent workflow nevertheless follows 
a relatively consistent series of steps: (i) fit Poisson nnodel, 
(ii) test other densities, (iii) develop structural nnodel on A, 
and (iv) refine interindividual variability nnodel. In a case of 



a parallel arnns placebo-controlled study, steps (i) and (ii) 
nnay be perfornned only on nontreated patients. Step (iii) 
includes the investigation of tinne effect, previous observation 
effect, covariate effect, drug effect, etc. Before step (iv), nnost 
paranneters nnay already consist in an association between 
a fixed and a randonn effect, so correlations (i.e., covariance 
nnatrices) and senniparannetric distributions (e.g., box-cox 
transformation) can be investigated. 

Model evaluation, as in any other nnodel exercise, takes 
place at each nnodel nnodification to quantify and qualify the 
fit innprovennent. 

An exannple in the fornn of a hands-on case study is pro- 
vided as Supplementary Material online. 

Model evaluation diagnostics 

Nunnerical evaluation in count nnodels is facilitated by the fact 
that nnost nnodels are nested with the Poisson nnodel. Hence, 
the likelihood ratio test can be applied between the hierarchical 
nnodels of a workflow following the steps previously described. 

Graphical evaluation does not include the sanne battery 
of goodness of fit plots as for continuous nnodels. Regard- 
less, residuals adapted to count nnodels are applicable.^^The 
Pearson residual, based on the estinnated nnean and variance 
connpared with the observations, and the Ansconnbe residual, 
built only with the estinnated nnean and the observations, nnay 
be inspected vs. A or or tinne to detect nnisspecifications. 

Sinnulation-based evaluation rennains nonetheless the gold- 
standard due to its ability to reliably depict several aspects — 
variability and structure — of the nnodel efficiently. Based on a 
preferably high nunnber of sinnulations potentially including 
uncertainty a visual predictive check can, for instance, be pro- 
duced. A visual predictive check typically displays percentiles of 
a dependent variable plotted against an independent variable 
both for observations and sinnulations, overlap of fitted data and 
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Figure 7 Visual predictive check for count data. Each panel 
represents the probability of observing a range of counts against 
time {P{Y = 0),P(1 <Y< 2),P(3 <Y < 4),P(Y > 5)). The solid line is 
the median of the observed data and the ribbon the 95% confidence 
interval around the median (dotted line) of 500 simulations from a 
count model fitted to count data. 

sinnulated data is the goal. In the case of count data, the ordinate 
can be the nunnber of events although its integer nature nnay 
alter the plot. A favored representation consists in plotting the 
proportions of counts within certain ranges. In Figure 7, these 
proportions are featured with regard to tinne. Other graphs to 
consider to connplete the assessnnent include: a stratification on 
treatnnent and other innportant categorical covariates, a switch 
to dose on the x-axis if enough dose levels and to innportant 
continuous covariates, an adjustment of the count proportions 
to proportions of transition values between consecutive counts, 
and a replacennent of the y-axis by the variance-to-nnean ratio. 
Sonne software, like PsN,^^ feature an autonnation tool generat- 
ing the visual predictive check with sinnulations, autobinning, 
and percentile calculation being done with one connnnand.^^ 

Count data simulation 

Whether intended for nnodel evaluation or nnodel inferences, 
sinnulations are an innportant part of count data analysis. In 
software where the distribution function is predefined (e.g., 
R), it can be used both for estinnation and for sinnulation, 
whereas in software where a user-defined function is neces- 
sary (e.g., NONMEM), the code should be adapted between 
estinnation and sinnulation. 

The sinnulation code needed for Poisson distributed data 
is once again borrowed fronn the proportional hazard nnodel 
(Eq. 15). It follows Algorithnn 1, according to which the nunn- 
ber of tinnes n that r needs to be drawn fronn a unifornn distri- 
bution for the sunn X/f=o'°9 (''/c)/-^/ to reach 1 corresponds 
to the sinnulated nunnber of events. 




50 100 
Average daily concentration (ng/ml) 

Figure 8 Dose selection plot. The exposure-response is depicted 
as change from baseline in A at day 28 against average daily 
concentration. The blue elements are the model-based mean and 
95% prediction interval generated from 500 simulations on a fine 
grid. The two-dimensional densities, built from patients observations, 
indicate the dose groups from the last conducted studies. A placebo 
group had been studied and a placebo effect observed. 



t = 0; n=0; 
while t < 1 do 

r-Z^(l); 

t = t + log(r)/- A, 
if t < 1 then 
I n=n+l 
end 
end 



Algoritlim 1 : Sinnulation of Poisson randonn variable 



p = 0; n = 0; 
r-ZY(l); 
while p < r do 
p = p + P{Yij 
if p < r then 
I n=n+l 
end 
end 



n) 



Piy. = 0) = e-^' <^ log (r) = -A, <^ 



log (0 



(15) 



Algorithm 2: Sinnulation of count nnodel randonn variable 

For pnnfs different fronn the Poisson nnodel, a nnore gen- 
eral sinnulation code is necessary. This code is described in 
Algorithnn 2 and is suitable for all nnodels presented in this 
tutorial. 



CPT: Pharmacometrics & Systems Pliarmacology 



Modeling and Simulation of Count Data 

Plan 



11 



Count model-based drug development 

Count models are inevitably part of decision-making within the 
drug development landscape, which is why solid knowledge 
about the approach is necessary A fundamental aspect is that 
count modeling is a flexible approach detached from any par- 
ticular therapeutic area. As more data are analyzed and more 
knowledge acquired, it is nevertheless imaginable and desir- 
able that platform models are built in specific diseases. As of 
now, predominance of application of mixed-effects count mod- 
els seems to be situated in epilepsy^^'^^ and oncology.^'^ Knowl- 
edge about the biological processes involved is invaluable and 
enables the development of mechanistic components into the 
models, but insight into the physiological variations is at times 
insufficient, justifying the choice of empirical functions. 

The impact of count modeling can take place at all phases 
of drug development, as shown by analyzed studies rang- 
ing from animal experiments^^ to patient trials.^ Decisions 
that can be informed include go-nogo or dose selection. The 
function A is where the drug effect, subject of the evaluation, 
is most likely to be implemented and a target response to 
be defined. Graphical^^ representations should be sought 
(e.g., Figure 8) and parameter uncertainty considered when 
relevant. Drug development strategic questions generally 
involve several aspects such as efficacy and safety, in which 
case Gupta etal.^ combined a count model and a categorical 
model to predict the therapeutic index of a new formulation. 
When the different aspects concern the same end point, like 
frequency and severity of a unique type of events, models 
can be combined on another level, but this has not been 
done for count models yet. 

A determinant factor in clinical development is trial effi- 
ciency. Underpowered or poorly designed studies lead to 
nonconclusive and expensive trials. Important decisions in 
study design can be supported by stochastic simulations 
and estimations from a previously developed model. Next to 
sample size and study duration, the length of the time interval 
during which events are counted may be investigated in the 
case of count studies. Optimal design, recognized to signifi- 
cantly impact the degree of conclusiveness of an experiment 
through parameter precision, has been shown to be able 
to handle discrete/count data.^^"^^ Hence, from a vantage 
point of a drug developer, count modeling and simulation 
can enable more precise, more certain, faster, and cheaper 
decisions. 

CONCLUSIONS 

Modeling and Simulation is a powerful tool to characterize, 
quantify, understand, predict, power, optimize, and rational- 
ize (pre)clinical trial data and studies. As in the case of any 
tool, its impact is maximized if it is shaped to fit its purpose. In 
the case of count outcome, tools in the form of pmfs are avail- 
able. They come with assumptions as well as extensions, 
which were the topic of this tutorial. The Poisson model, a 
close relative of the survival model, is the basis for all count 
data models. Adaptations for handling overdispersion, under- 
dispersion, autocorrelation, or inhomogeneity were pro- 
posed in the literature and presented here. Other sources of 
information may be useful to the reader, such as Cameron 



etal.^^ or Coxe etal..^^ All in all, count modeling is an integral 
part of the dynamic science of pharmacometrics alongside 
count data being possible end points collected during drug 
development. 

Acknowledgments. We thank Mats Karlsson for encourag- 
ing data type-informed modeling analyses and Sebastian 
Ueckert for valuable comments. 

Conflict of Interest. The author declared no conflict of interest. 

1 . Barrett, J.S., Fossler, M.J., Cadieu, K.D. & Gastonguay, M.R. Pharmacometrics: a 
multidisciplinary field to facilitate critical thinking in drug development and translational 
research settings. J. Clin. Pharmacol. 48, 632-649 (2008). 

2. Rowland, M., Benet, L.Z. & Riegelman, S. Two-compartment model for a drug and its 
metabolite: application to acetylsalicylic acid pharmacokinetics. J. Pharm. Sci. 59, 364-367 
(1970). 

3. Sheiner, L.B. & Beal, S.L. Evaluation of methods for estimating population 
pharmacokinetics parameters. I. Michaelis-Menten model: routine clinical pharmacokinetic 
data. J. Pharmacoklnet. Blopharm. 8, 553-571 (1980). 

4. Paule, I., Girard, R, Freyer, G. & Tod, M. Pharmacodynamic models for discrete data. Clin. 
Pharmacoklnet. 51 , 767-786 (201 2). 

5. Mould, D.R. & Upton, R.N. Basic concepts in population modeling, simulation, and model- 
based drug development. CPT Pharmacometrics Syst. Pharmacol. ^, e6 (2012). 

6. Holford, N. A time to event tutorial for pharmacometricians. CPT: Pharmacometrics Syst. 
Pharmacol!, e43 (2013). 

7. Frame, B., Miller, R. & Lalonde, R.L. Evaluation of mixture modeling with count data using 
NONMEM. J. Pharmacoklnet. Pharmacodynamics 30, 167-183 (2003). 

8. Gupta, S.K., Sathyan, G., Lindemulder, E.A., Ho, PL, Sheiner, L.B. & Aarons, L. 
Quantitative characterization of therapeutic index: application of mixed-effects modeling 
to evaluate oxybutynin dose-efficacy and dose-side effect relationships. Clin. Pharmacol. 
Ther65, 672-684 (1999). 

9. Wirz, S., Wartenberg, H.C. & Nadstawek, J. Less nausea, emesis, and constipation 
comparing hydromorphone and morphine? A prospective open-labeled investigation on 
cancer pain. Support. Care Cancer ^6, 999-1009 (2008). 

1 0. Sheiner, L.B. The population approach to pharmacokinetic data analysis: rationale and 
standard data analysis methods. Drug Metabolism Rev. 15, 153-171 (1984). 

1 1 . Plan, E.L., Ma, G., Nagard, M., Jensen, J. & Karlsson, M.O. Transient lower 
esophageal sphincter relaxation pharmacokinetic-pharmacodynamic modeling: 
count model and repeated time-to-event model. J. Pharmacol. Experiment. Then 339, 
878-885 (2011). 

12. Simeon Denis Poisson and Christian Heinrich Schnuse. Recherches sur la Probabilite des 
Jugements en matiere criminelle et en matiere civile. (Meyer, Parris, France, 1841). 

1 3. de Moivre, A. De mensura sorits seu, de probabilitate eventuum in ludis a casu fortuito 
pendentibus. (1711). 

14. Holford, TR.The analysis of rates and of survivorship using log-linear models. Biometrics 
36, 299-305(1980). 

1 5. Laird, N. & Olivier, D. Covariance analysis of censored survival data using log-linear 
analysis techniques. J. Am. Stat. Association 76, 231-240 (1981). 

1 6. Upton, R.N. & Mould, D.R. Basic concepts in population modeling, simulation, and model- 
based drug development: part 3-introduction to pharmacodynamic modeling methods. 
CPT Pharmacometrics Syst. Pharmacol. 3, e88 (2014). 

1 7. Mould, D.R. & Upton, R.N. Basic concepts in population modeling, simulation, and model- 
based drug development-part 2: introduction to pharmacokinetic modeling methods. CPT 
Pharmacometrics Syst. Pharmacol 2, e38 (2013). 

18. Hoel, PG. On indices of dispersion. Ann. Math. Stat. 14, 155-162 (1943). 

1 9. Selby B. The index of dispersion as a test statistic. Blometrlka 52, 627-629 (1 965). 

20. Hoaglin, D.C. A Poissonness plot. The American Statistician 34, 146-149 (1980). 

21 . Hoaglin, D.C, Mosteller, F & Tukey J.W. Exploring Data Tables, Trends, and Shapes. (John 
Wiley & Sons, Hoboken, NJ, 2011). 

22. Gardner, W., Mulvey E.P & Shaw, E.G. Regression analyses of counts and rates: Poisson, 
overdispersed Poisson, and negative binomial models. Psychol. 118, 392-404 
(1995). 

23. Land, K.C., McCall, PL. & Nagin, D.S. A comparison of Poisson, negative binomial, and 
semiparametric mixed Poisson regression models with empirical applications to criminal 
careers data. Sociological Methods & Research 24, 387^42 (1996). 

24. Troconiz, I.F, Plan, E.L., Miller, R. & Karlsson, M.O. Modelling overdispersion and 
Markovian features in count data. J. Pharmacoklnet. Pharmacodyn. 36, 461-477 
(2009). 

25. Ahn, J.E., Plan, E.L, Karlsson, M.O. & Miller, R. Modeling longitudinal daily seizure 
frequency data from pregabalin add-on treatment. J. Clin. Pharmacol. 52, 880-892 
(2012). 

26. Lloyd-Smith, J.O., Schreiber, S.J., Kopp, PE. & Getz, W.M. Superspreading and the effect 
of individual variation on disease emergence. Nature 438, 355-359 (2005). 



www.nature.com/psp 



Modeling and Simulation of Count Data 

Plan 



12 



27. Consul, RC. & Jain, G.C. A generalization of the Poisson distribution. Technometrics 1 5, 
791 (1973). 

28. Consul, PC. & Shoukri, M.M. The generalized poisson distribution when the sample 
mean is larger than the sample variance. Communications in Statistics - Simulation and 
Computation U, 667-681 (1985). 

29. Efron, B. Double exponential families and their use in generalized linear regression. J. Am. 
Stat. Association 81 , 709 (1 986). 

30. Ridout, M.S. & Besbeas, R An empirical model for underdispersed count data. Statistical 
Modeiiing^, 77-89 (2004). 

31 . Bonate, PL., Sung, C, Welch, K. & Richards, S. Conditional modeling of antibody 
titers using a zero-inflated poisson random effects model: application to Fabrazyme. J. 
Pliarmacol<inet Pliarmacodyn. 36, 443-459 (2009). 

32. Ambagaspitiya, R.S. & Balakrishnan, N. On the compound generalized Poisson 
distributions. Astin Buiietin 24, 255-263 (1994). 

33. King, G. Variance specification in event count models: from restrictive assumptions to a 
generalized estimator. Am. J. Poiiticai ScL 33, 762-784 (1 989). 

34. Winkelmann, R. & Zimmermann, K.F. Count data models for demographic data. Matli. 
Popui. Stud 4, 205-221 , 223 (1 994). 

35. Zou, Y, Lord, D. & Geedipally S.R. Over- and Under-Dispersed Count Data: Comparing 
Conway-Maxwell-Poisson and Double-Poisson Distributions. Transportation Research 
Board 91st Annual Meeting, Washington, DC, USA, 2012. 

36. Zou, Y, Geedipally S.R. & Lord, D. Evaluating the Double Poisson Generalized Linear 
Model. Transportation Research Board 92nd Annual Meeting, Washington, DC, USA, 2013. 

37. Plan, E.L, Elshoff, J.P, Stockis, A., Sargentini-Maier, M.L & Karlsson, M.O. Likert pain 
score modeling: a Markov integer model and an autoregressive continuous model, din. 
Pliarmacoi. Tlier. 91 , 820-828 (2012). 

38. Karlsson, M.O., Beal, S.L. & Sheiner, L.B. Three new residual error models for population 
PK/PD analyses. J. Pliarmacoldnet. Biopliarmaceutics23, 651-672 (1995). 

39. Tornoe, C.W., Overgaard, R.V., Agerso, H., Nielsen, H.A., Madsen, H. & Jonsson, E.N. 
Stochastic differential equations in NONMEM®: implementation, application, and 
comparison with ordinary differential equations. Pliarmaceuticai Res. 22, 1247-1258 
(2005). 

40. Baum, L.E. & Petrie, T Statistical inference for probabilistic functions of finite state Markov 
chains. Ann. Math. Stat. 37, 1 554-1 563 (1 966). 

41 . Delattre, M., Savic, R.M., Miller, R., Karlsson, M.O. & Lavielle, M. Analysis of exposure- 
response of CI-945 in patients with epilepsy: application of novel mixed hidden Markov 
modeling methodology J. Pliarmacol<inet. Pliarmacodyn. 39, 263-271 (2012). 

42. Zeger, S.L. & Qaqish, B. Markov regression models for time series: a quasi-likelihood 
approach. Biometrics 1019 (1988). 

43. Velez de Mendizabal, N., Troconiz, I.R, Hutmacher, M.M. & Bies, R. Discrete distribution 
models for relapsing-remitting dynamics observed in multiple sclerosis. PAGE2^ (2012). 

44. Schoemaker, R., Snoeck, E. & Stockis, A. Exposure-response modeling of daily seizure 
counts in focal epilepsy trials. PAGE 19 (2010). 

45. Singpurwalla, N.D. The hazard potential. J. Am. Stat. Association 1 01 , 1 705-1 71 7 (2006). 

46. Cox, E.H., Veyrat-Follet, C, Beal, S.L., Fuseau, E., Kenkare, S. & Sheiner, L.B. A 
population pharmacokinetic-pharmacodynamic analysis of repeated measures 
time-to-event pharmacodynamic responses: the antiemetic effect of ondansetron. J. 
Pharmacokinet. Biopharmaceutics 27, 625-644 (1999). 

47. Simeoni, M. et al. Counting events: population approaches using NONMEM. PAGE 14, 
2005. 

48. Plan, E.L, Karlsson, K.E. & Karlsson, M.O. Approaches to simultaneous analysis of 
frequency and severity of symptoms. Clin. Pharmacol. Ther 88, 255-259 (2010). 



49. Yano, I., Beal, S.L. & Sheiner, L.B. The need for mixed-effects modeling with population 
dichotomous data. J. Pharmacokinet. Pharmacodyn. 28, 389^12 (2001). 

50. Plan, E.L., Maloney A., Troconiz, I.F & Karlsson, M.O. Performance in population 
models for count data, part I: maximum likelihood approximations. J. Pharmacokinet. 
Pharmacodyn. 36, 353-366 (2009). 

51 . Savic, R. & Lavielle, M. Performance in population models for count data, part II: a new 
SAEM algorithm. J. Pharmacokinet. Pharmacodyn. 36, 367-379 (2009). 

52. Metrum Institute. MI255: Exposure-Response Modeling of Categorical, Count, and Time- 
to-Event Data Using Bayesian Methods, <http://www.youtube.com/metruminst>. (2011). 
Lecture 6. 

53. Pierce, D.A. & Schafer, D.W. Residuals in generalized linear models. J. Am. Statistical 
Association 81, 977-986 (1986). 

54. Lindbom, L, Pihigren, P, Jonsson, E.N. & Jonsson, N. PsN-Toolkit-a collection of 
computer intensive statistical methods for non-linear mixed effect modeling using 
NONMEM. Computer Methods and Programs in Biomedicine 79, 241-257 (2005). 

55. Keizer, R.J., Karlsson, M.O. & Hooker, A. Modeling and Simulation Workbench for 
NONMEM: Tutorial on Pirana, PsN, and Xpose. CPT Pharmacometrics Syst. Pharmacol 2, 
e50 (2013). 

56. Jonker, D.M., van de Mheen, C, Filers, PH., Kruk, M.R., Voskuyl, R.A. & Danhof, M. 
Anticonvulsant drugs differentially suppress individual ictal signs: a pharmacokinetic/ 
pharmacodynamic analysis in the cortical stimulation model in the rat. Behav. NeuroscL 
117,1076-1085 (2003). 

57. Ito, K. & Murphy D. Application of ggplot2 to Pharmacometric Graphics. CPT 
Pharmacometrics Syst. Pharmacol 2, e79 (201 3). 

58. Box, J.F R.A. Fisher and the design of experiments, 1 922-1 926. The American Statistician 
34,1-7(1980). 

59. Nestorov, I., Graham, G., Duffull, S., Aarons, L., Fuseau, E. & Coates, P Modeling and 
stimulation for clinical trial design involving a categorical response: a phase II case study 
with naratriptan. Pharm. Res. 18, 1210-1219 (2001). 

60. Ogungbenro, K. & Aarons, L. Sample size calculations for population pharmacodynamic 
experiments involving repeated dichotomous observations. J. Biopharmaceutical Stat. 18, 
1212-1227 (2008). 

61 . Nyberg, J., Karlsson, M.O. & Hooker, A. Population optimal experimental design for 
discrete type data. PAGE^S (2009). 

62. Cameron, A.C. & Trivedi, PK. Regression Analysis of Count Data. (Cambridge University 
Press, Cambridge, England, 1998). 

63. Coxe, S., West, S.G. & Aiken, L.S. The analysis of count data: a gentle introduction to 
poisson regression and its alternatives. J. Personality Assessment 9^ , 121-136 (2009). 



Commons 
Unpolled 



rJci®®© ^^'^ ^^^^ licensed under a Creative 
Attribution-NonCommercial-NoDerivs 3.0 



License. The images or other third party material in this article are 
included in the article's Creative Commons license, unless indicated 
otherwise in the credit line; if the material is not included under the 
Creative Commons license, users will need to obtain permission from 
the license holder to reproduce the material. To view a copy of this 
license, visit http://creativecommons.orc|/licenses/by-nc-nd/3.0/ 



Supplementary information accompanies this paper on the CPT: Pharmacometrics & Systems Pharmacology \Nebs\{e 

(http://www.nature.com/psp) 



CPT: Pharmacometrics & Systems Pharmacology 



